#include "base.h"

int main(int argc,char **argv)
{

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

          //////////////////
          // Input   ///////
          //////////////////

  char dir[1024];

   char setpprint[16];
   char Oa2print[16];
   char output_name[144];
   int list_prop[16];
   int Nprop;
   int beta100;
   int Nmass, imass; 
   char momentum_order[1024];

          ////// Parsing command line


          for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
            if(strncmp(argv[i],"-dir",4)==0)  sscanf(argv[i],"-dir=%s",dir);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
            if(strncmp(argv[i],"-beta100",8)==0)  sscanf(argv[i],"-beta100=%d",&beta100);
            if(strncmp(argv[i],"-Nmass",6)==0)  sscanf(argv[i],"-Nmass=%d",&Nmass);

            if(strncmp(argv[i],"-momentum_order",15)==0)  sscanf(argv[i],"-momentum_order=%s",momentum_order);
	    if(strncmp(argv[i],"-Nprop",6)==0)  sscanf(argv[i],"-Nprop=%d",&Nprop);
	    if(strncmp(argv[i],"-output_name",6)==0)  sscanf(argv[i],"-output_name=%s",output_name);
            if(strncmp(argv[i],"-list_prop",10)==0)  sscanf(argv[i],"-list_prop=%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d",&list_prop[0],&list_prop[1],&list_prop[2],&list_prop[3],&list_prop[4],&list_prop[5],&list_prop[6],&list_prop[7],&list_prop[8],&list_prop[9],&list_prop[10],&list_prop[11],&list_prop[12],&list_prop[13],&list_prop[14],&list_prop[15]);
            
	    if(strncmp(argv[i],"-Oa2",4)==0)  sscanf(argv[i],"-Oa2=%s",Oa2print);


          }
   int Nmass_comb=Nmass;
   for (int i=1; i<Nmass; i++)  Nmass_comb=Nmass_comb+(Nmass-i);

   // 
   printf("Promedium over %d democratic sets: \n",Nprop);
  for (int iprop=0; iprop<Nprop; iprop++)
  {
    whichProp(list_prop[iprop], setpprint);
    printf("%s ",setpprint);
  }
  printf("\n%s\n",output_name);


   /////////////////////////
  //// Momenta set
  //////////////////////////////
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom_%d.dat",beta100);
  read_momenta_data(momentafile);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(momentum_order,nn);
  // Remember that all set_p corresponding to each iprop must have the same p2. This is true if I use the even (xxx+) or the odd (xxx-) set separetely
  Momentum(momentum_order);
  // equivalent momenta
  Nmom_eq=Number_of_Equivalent_Momenta(SinP2);
  Average_over_Equivalent_Momenta(SinP2,SinP2,SinP2_eq);

	// IO

  char file[1024];
  char z_output[1024];
  char dir_qav[1024];
  char gammaprint[8];

       //////////////////////////////////
        // Read Zq and do promedium
        /////////////////////////////



         double ***Zq=malloc(sizeof(double**)*(Nclust+1));
         for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust]=malloc(sizeof(double*)*(Nmom_eq));
	 for (int iclust=0; iclust<Nclust+1; iclust++)  for (int imom=0; imom<Nmom_eq; imom++) Zq[iclust][imom]=malloc(sizeof(double)*(Nmass));

         double ****ZO=malloc(sizeof(double***)*(Nclust+1));
         for (int iclust=0; iclust<Nclust+1; iclust++) ZO[iclust]=malloc(sizeof(double**)*(Nmom_eq));
	 for (int iclust=0; iclust<Nclust+1; iclust++)  for (int imom=0; imom<Nmom_eq; imom++) ZO[iclust][imom]=malloc(sizeof(double*)*(Nmass_comb));
         for (int iclust=0; iclust<Nclust+1; iclust++)  for (int imom=0; imom<Nmom_eq; imom++) for (int imass=0; imass<Nmass_comb; imass++) ZO[iclust][imom][imass]=malloc(sizeof(double)*(5));

        double **aux_allmomenta=malloc(sizeof(double*)*(Nmom_eq));
       for (int imom=0; imom<Nmom_eq; imom++) aux_allmomenta[imom]=malloc(sizeof(double)*(Nclust+1));

        double *aux_clust=malloc(sizeof(double)*(Nclust+1));

 // 1. Initialize to zero
  for (int iclust=0; iclust<Nclust+1; iclust++) for (int imom=0; imom<Nmom_eq; imom++) for (int imass=0; imass<Nmass; imass++)
  {
	Zq[iclust][imom][imass]=0;
   }
  for (int iclust=0; iclust<Nclust+1; iclust++) for (int imom=0; imom<Nmom_eq; imom++) for (int imass=0; imass<Nmass_comb; imass++)
  {
	for (int igamma=0; igamma<5; igamma++) ZO[iclust][imom][imass][igamma]=0;
  }
 // 2. Read
	printf("Reading \n");
  for (int iprop=0; iprop<Nprop; iprop++)
  {


    whichProp(list_prop[iprop], setpprint); 
    sprintf(dir_qav,"%s/%s/QAV/",dir,setpprint);
	printf("%s\n",setpprint);
    // Zq
    for (int imass=0; imass<Nmass; imass++) {
    sprintf(file,"Zq_%s_m%d",Oa2print,imass);
    for (int imom=0; imom<Nmom_eq; imom++)
    {
           read_input_onemomenta(dir_qav,file,aux_clust,imom);
           for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust][imom][imass]=Zq[iclust][imom][imass]+aux_clust[iclust];
    }
    }

   
    // ZO 
   char gammaprint[8];
   imass=0;
   for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
   {
   for (int igamma=0; igamma<5; igamma++)
   {
         whichGammaProm(igamma,gammaprint);
	 sprintf(file,"Z%s_%s_m%d%d",gammaprint,Oa2print,im1,im2);
            for (int imom=0; imom<Nmom_eq; imom++)
            {
                read_input_onemomenta(dir_qav,file,aux_clust,imom);
                for (int iclust=0; iclust<Nclust+1; iclust++) ZO[iclust][imom][imass][igamma]=ZO[iclust][imom][imass][igamma]+aux_clust[iclust];
            }

    }
   imass++;
   }
   }
  // 3. Divide by the number of democratic sets considered
  for (int iclust=0; iclust<Nclust+1; iclust++) for (int imom=0; imom<Nmom_eq; imom++) for (int imass=0; imass<Nmass; imass++)
  {
        Zq[iclust][imom][imass]=Zq[iclust][imom][imass]/Nprop;
  }
  for (int iclust=0; iclust<Nclust+1; iclust++) for (int imom=0; imom<Nmom_eq; imom++) for (int imass=0; imass<Nmass_comb; imass++)
  {
        for (int igamma=0; igamma<5; igamma++) ZO[iclust][imom][imass][igamma]=ZO[iclust][imom][imass][igamma]/Nprop;
  }




        /////////////////
        // Write output
        ////////////////

    sprintf(dir_qav,"%s/%s/QAV/",dir,output_name);
    printf(" Writting output in  %s\n",dir_qav);

	double **aux_allmomenta_real=malloc(sizeof(double*)*(Nmom_eq));
	for (int imom=0; imom<Nmom_eq; imom++) aux_allmomenta_real[imom]=malloc(sizeof(double)*(Nclust+1));

      for (int imass=0; imass<Nmass; imass++)
      {
         sprintf(z_output,"%s/Z_ALLP/Zq_%s_m%d.dat",dir_qav,Oa2print,imass);
         for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=Zq[iclust][imom][imass];
         write_output_allmomenta_real(z_output,aux_allmomenta_real,SinP2_eq,Nmom_eq);

         sprintf(file,"Zq_%s_m%d",Oa2print,imass);
         write_output_momentabymomenta_real(dir_qav,file,aux_allmomenta_real,SinP2_eq,Nmom_eq);
       }

   imass=0;
   for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
   {
   for (int igamma=0; igamma<5; igamma++)
   {
        whichGammaProm(igamma,gammaprint);

         sprintf(z_output,"%s/Z_ALLP/Z%s_%s_m%d%d.dat",dir_qav,gammaprint,Oa2print,im1,im2);
         for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=ZO[iclust][imom][imass][igamma];
         write_output_allmomenta_real(z_output,aux_allmomenta_real,SinP2_eq,Nmom_eq);


          sprintf(file,"Z%s_%s_m%d%d",gammaprint,Oa2print,im1,im2);
         write_output_momentabymomenta_real(dir_qav,file,aux_allmomenta_real,SinP2_eq,Nmom_eq);

   }
   imass++;
   }

	free(Zq); free(ZO);
    return 0;

}

